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ABSTRACT 


The objective of this thesis is the refinement of a method for 
evaluating the pressure distribution on a body surface of arbitrary 
shape in incompressible flow. 

The solution to the problem is obtained in terms of the velocity 
potential. However, the formulation for the potential cannot be 
solved in closed form with the exception of a few special cases. 

Hence, the solution must be obtained through numerical approximations 
which require the use of a high speed digital computer. 

Two different numerical formulations are considered here. In 
the first method, the box method , the body surface is approximated by 
N small elements, 2 — . Within these elements or boxes, the value of 
the potential, <p , is considered constant and equal to the value at 
the center of the element, . This results in a system of N linear, 
algebraic equations in N unknowns. 

The other method is the collocation method in which the unknown 
function is approximated by a linear combination of N prescribed 
functions which satisfy the boundary conditions of the problem and 
N unknown coefficients resulting in a system of linear, algebraic 
equations. The specific approach applied here is called the modal 
method . 

In this thesis, both of these methods are applied to a very thin, 
rectangular wing in incompressible, steady flow. Of the two methods. 
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the box method is more practical as it is applicable to more general 
geometries (the modal method requires a new set of functions for each 
geometry) and requires less computer time (fifty percent of that 
required by the modal method for the same problem). 
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1. THEORETICAL FORMULATION 


1.1 Introduction 

The problem of evaluating the pressure distribution on a wing 
immersed in a fluid has, in the past, and is presently attracting a 
great deal of attention. The most conmon method applied is that of 
lifting surface theories using the velocity potential. The principal 
equation used is an integral equation which relates the pressure- 
difference distribution over the body surface with the downwash pro- 
duced by the body surface. This integral equation may be expressed as^. 

where \.he pressure difference, » is the unknown, and 
Is the so-called Kernel function. 

There are two major objections to this type of approach. The first 
of these is that the numerical solution of this equation is quite 
complicated. The difficulties arise from the rather complicated Kernel 
function and the numerical integration of improper integrals. The second 
objection is that lifting-surface theories are not easily generalized to 
Include more complicated geometries. Similar objections hold for the 
problem of oscillating wings. An excellent analysis of the status of art 
is given in Ref. 2. 
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In order to circumvent these problems, a general method has been 
developed in Ref. 3. This method is valid for both compressible and 
incompressible potential flows around bodies of arbitrary shape in 
arbitrary motion. The formulation is based on Green's formula for the 
equation of unsteady compressible potential aerodynamic flow. 

Although theoretical formulation given in Ref. 3 is very general , the 
application of the method is limited to thin wings in steady subsonic 
flow. The objective of this thesis is to explore further the applica- 
bility of the method mentioned above. However, this analysis is 
limited, for simplicity, to steady, incompressible flow. In the rest 
of this section, the method developed in Ref. 3 is summarized for the 
particular case of steady, incompressible flow. 

1.2 Statement of the Problem 

The equation for steady incompressible potential aerodynamic flow 
is given by the well known Laplace's equation^ 


= ° 


1.2 


where 0 is the velocity potential. Let the surface of the body, 
be given by the equation 

s B (W,z)^ o 
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t 



Then the tangency boundary conditions are given by 


$0 

&Y\ 



1.4 


where ri is the normal to the surface 21. Finally the boundary condi- 
tion at co is given by 


0 - UooX 


1.5 


By introducing the perturbation potential, zp » such that 


0 = U* (X + <t>) 

equations 1.2, 1.4, and 1.5 yield Eqs. 1.7, 1.8, and 1.9 respectively 
which are given as: 



1.7 


and 


_ 1 b S* 

OO b X 


on the surface 



and finally 


1.8 


9=0 


1.9 


at Infinity. The problem is to determine the solution cf Eq. 1.7 using 
the boundary conditions given by Eqs. 1.8 and 1.9. Finally the pressure 
coefficient may be obtained from the linearized Bernoulli Theorem as^ 


P-P<*> = 



2 . 

CD cTx 


1.10 


and the pressure coefficient is given by 



P Poo 



i.ii 


1.3 Formulation of the Problem 

The basis for the formulation of the problem is Green's formula for 
the Laplace equation. Green's formula can be expressed as^ 


Q (X A X) =~jf ■$*; (*&!,) ^ <1X 

St r J d X 


1.12 
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with 


1 

=s Q*-x 1 f+ { v -y,f + -2-^3 2 


where (x, y, z) is the control point and (x 1 , y^ , z^) is the dummy 
point of integration on the surface X* The surface surrounds the 
body and the wake. 

From Eq. 1.13 it can be seen that the potential is given by an 
integral representation composed of sources , t / A 'It' P , of intensity 
(£><?/ dti ) and of doublets, b/d\\ (1/j-Tc' IT ) , of intensity <Q . The 
values of are known from the boundary conditions as shown by 

Fig. 1. All that remains to be determined is the value of Cp which is 
unknown on the body surface. 

In order to evaluate <p on the surface X. let the control point 
approach a point on the surface ]T . As shown in Appendix A this 
yields 




1.14 


+ # <SY, ,2,)^ (—Hell 
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Note that Eq. 1.14 differs from Eq. 1.13 by the factor 1/27 f Instead of 
1/4 TT . In addition the control point (x, y, z) is on the surface 
instead of outside. If is known everywhere on 3L » Eq. 1.14 is a 

Freedholm integral equation of second kind, whose solution exists and is 
unique (exterior Von Newman's problem of Laplace's equation 0 ). However, 
the surface surrounds not only the body but also the wake. The procedure 
for treating the contribution of the wake is discussed in the next Sub- 
section, where the integral equation given by Eq. 1.15 is obtained and 
Che methods of solution are outlined in Subsection 1.5. 

1.4 The Vfake Contribution 

As mentioned above, the surface X. in Eq. 1.14 is composed of two 
branches; the first,lg» surrounds the body and the second,X w , surrounds 
the wake. The treatment of the wake introduced in Ref. 3 is used here. 

The wake is assumed to be composed of straight vortex-lines emanating from 
the trailing edge of the body. Under, this assumption. Eq. 1.14 becomes 
(see Appendix B) 



with the wake integral, I w , given by 


1.15 
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1.16 


Iw=jj Jr(wr) dI 

where IT w is the upper side of the surface of the wake. I.i the problem 
considered in the following, the trailing edge is on the plane Z = 0 
and thus Eq. 1.16 reduces to 

T w — \ 

where b is the span of the wing and J is given by 

W 



Note that the wake formulation is based on bodies having sharp trailing 
edge otherwise it is not possible to know where the wake leaves the 
body. 
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1.5 Methods. of Solution 


The necessary relation for evaluation cp on the body surface 
is given in Subsection 1.4 as shown in Eq. 1.15. It may be noted that, 
with the exception of a few special cases, it cannot be solved in 
closed form. Hence the solution must be obtained through numerical 
approximate methods, which generally require the use of a high speed 
digital computer. 

Several methods can be used. The one used in Ref. 3, which will 
be called box method , is summarized here. The surface is divided 

into N small elements (boxes) £ • . The value of is assumed to be 

constant within each element and equal to the value ^ at the center 

of the element. By satisfying Eq. 1.15 at the centers of the elements 
^ , one obtains a system of N linear algebraic equations with N 
unknowns, , which can be easily solved on a computer. A more 
detailed description of the box method is given in Section 2. 

However, one of the most successful methods for lifting surface 
theories is the collocation method in which the unknown function is 
approximated by a linear combination of N prescribed functions with N 
unknown coefficients. The prescribed functions satisfy the boundary 
conditions of the problem whose N unknown coefficients are determined 
by solving at N given points (collocation points). Thus before 
proceeding into a more elaborated problem, the legitimate question of 
which of the two methods is more convenient for this problem must be 
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answered. A method equivalent to the collocation method (calle d modal 
method ) is described in Section 3. The advantage to this method is chat 
it gives a smooth curve solution rather than a solution at discrete 
points. However, it is limited in that it requires a different set of 
functions for each particular geometry with increasing difficulty in 
guessing appropriate functions with more complex geometries, e.g., wing- 
body interference. Also, this method requires large amounts of computer 
time, nearly twice the amount required by the other methods explored. 

Thus it was decided to continue the analysis of the problem by 
using the box method. Next it should be noted that in Ref. 3 the 
coefficient of the equations of the box method are evaluated analytically 
by assuming that each surface element can be rtplaced by its tangent 
plane at the center of the element. This hypothesis is adequate for very 
thin wings (thickness ratio 'Y £1%). However, even for moderately thick 
wings ('Y =10%) the condition that, in each equation, the sum of the 
coefficients be equal to one is only poorly satisfied’ . Thus another 
legitimate question must be answered: can the convergence of the solution 

be improved by using the correct surface element and evaluating the 
coefficients numerically (instead of analytically)? This method is 
described in detail in Section 4. The results obtained in Section 6, 
however, indicate that the condition that in each equation the sum of 
the coefficients be equal one is fairly well satisfied for thick wings, 
but poorly satisfied for thin wings. 
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Hence, a third approach for the evaluation of the coefficients 
was explored. The coefficients are first evaluated analytically with 
the tangent plane hypothesis. Then they are corrected by adding the 
Integral of the difference (between correct-element integrand and 
tangent-plane integrand) evaluated numerically. The details of this 
procedure are given in Section 5. The results (Section 6) indicate 
that this procedure (numerical-analytic evaluation of coefficients) is 
better than the purely analytic and the purely numerical procedure. 

A diagram outlining the development of the above methods is given in 
Fig. 2. 
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2. THE BOX METHOD 


2.1 Geometry of the Problem 

The body used for studying the relative efficiency of the 
different methods is the same as that used in Ref. 3, i.e., a rec- 
tangular* thin wing in steady flow. For this particular problem, 
alternate numerical solutions (lifting surface theory. Ref. 4) as well 
as experimental results (Ref. 5) are available. 

The wing used has a symmetric profile with thickness h given by 


h = 



- *\ 2.1 


The thickness ratio (i.e., the ratio of the maximum body thickness to 
the chord) is given by 




2.2 


Rearranging Eq. 2.1 gives the geometry of the wing at zero angle of 
attack in terms of ^ and vj as 
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X - c I 

y - M 

z = t rc b 

with ^ and ^ given by 



and plus (minus) is for the upper (lower) side in Eq. 2.6. 

For angle of attack, , different from zero, the following 
coordinate transformation is used 


— ■ 

- $ Hi -n 


2.3 


X =• K CoS + 2: Stn<K 

y = 7 

'Z = - X Sm®< + 2 
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For purposes of comparison with the methods presented here and other 
methods, all results are evaluated at an angle of attack of 5° (Fig. 3). 
2.2 The Box Method 

The box method is used to evaluate the coefficients of the source, 
doublet and wake contributions. First consider the body surface, £ 
as divided into elements d£.j. this assumption, Eq. 1.15 can 

be written as 



Applying the Mean Value Theorem and assuming that the potential can be 
approximated by the value at the center of the element yields 
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where the index applies to the trailing edge boxes only for wake 
considerations. Fi-'lly, satisfying Eq. 2.8 at P(k), the centers of 
the elements .yield 


^ k. x Qi -+-1E. q.i 


The coefficients b^, c^, and w^ are defined by 



2.10 



Rearrangment gives the final form as 





= O 


2.11 


with 


a 


VU 




2.12 


where S ^ is the usual Kronecker Delta. 

2.3 Definition of Boxes 

In defining the boxes on the body surface, first consider the 
classical transformations. 


x - c 


f ^ 1 
-i - +1 


2.13 


y - 


-in 


This reduces the wing to one with chord equal to one and span equal to 
two (Fig. 4). iiote that ^3 / , and £5 /c)y^ are infinite at the 

leadiny edge and tip of the wing, respectively. This suggests a modifi- 
cation of the body surface boxes that gives a better description of the 
leading edge and the tip. This is best accomplished by using a non- 
uniform mesh for the definition of the boxes which is expressed mathe- 
matically by the transformation. 



O £ X £ 1 

-1 £ 9 £ i 


2.14 


This causes the boxes at the leading edge and the tip to be smaller 
than those not along the leading edge and 'tin -Fig. 5). Boxes of 
constant size AX, AY in the plane X, V given by 
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AX - 1/ N^X 

2.15 

Ay = 1 /hiy 

are used where NX and MY are the number of boxes in direction X and Y, 
respectively. In other words, the center of the box, (m, n) is given by 
(see Fig. 5), 

_(C) 

2.16 

(C) -< — 

V = I'n- })*v w- 1,... ,wy 

whereas, its boundaries are given by. 


- <*0 


_ <p> 

v £ 

x 

- 



-CP) 

y - 

y 

± y 
7 r \ 


with 


l 


— (tn) 

- cel 


Ay _ 

C m - 


X = 

vn 

X 

tr> 


a 

1 } A X 


-(C) 

-4- 

AX 



_ 

tn 


z 

- rr\ 

AX 

— Itn 

- Cc) 

— 

II 

( n- 

1 )^y 

— tp> 

Yt* " 

^*1 


II 

in a9 


Note that, for each set of values m and n, there are two boxes, one 
on the upper and one on the lower side of the wing; hence, the total 
number of boxes on the right hand part of the wing is 

n - a. mx * ny 2.19 

2.4 Comments 

The results obtained in Ref. 3 indicate that the method is both 
fast and accurate for very small thicknesses. However, the condition 
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(which should bo satisfied by the coefficients Cj,^ , Ref. 3, Eq. 7.22) is 
satisfied within a 4% tolerance for thickness ratio r = 1% while for 
T= 10%, the condition i only satisfied within a 24% tolerance. In 
order to improve this condition, two approaches are available: trying 
a different method or improving the evaluation of the coefficients in 
the box method. Both approaches are explores in Section 3 and Section 4 
respectively. 



.3. THE MODAL METHOD 


3.1 Introduction 

The so-called modal method is a method for determining the 
pressure distribution on a body surface by approximating the shape of 
the aerodynamic potential with a series of mode shapes in the spanwise 
and chordwise directions. Once the shape of the potential has been 
determined, the pressure distribution can be obtained by differentiating 
according to Eq. 1.11. The results obtained from this method are in 
very close agreement with those obtained with the box method. 

Although a straight forward, fairly uncomplicated method that 
gives good results in a wide range of thickness ratios, this method 
has one undesirable quality. It is quite good for relatively simple geome- 
tries but assuming the series representation of the potential for more 
involved geometries is quite difficult. One of the main interests in 
evaluating the pressure distribution is for wing-body interference and 
for this method, estimating a potential shape such as this would be 
quite difficult. Also, this method requires large computer time, up to 
twice as much as the box method. 

3.2 Theory 

As mentioned in Section 1, an alternate method, the modal method, 
is explored before attacking the problem of thick wings. The modal method 
approach is to approximate the shape of the potential by a series of mode 
shapes in the chordwise and spanwise directions. This results in a 
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system of N mode equations with N unknown coefficients. Once the 
potential distribution has been determined, the pressure distribution 
can be determined using Eq. 1.11. 


The modes of cp are chosen as follows: The function <£> can be 
separated into symmetric and anti-symmetric components with respect to 
the midplane of the wing (Fig. 6a & b) as. 


<P 



b <£ 


A 


3.1 


The symmetric part is given by (Fig. 6a), 



P Q 


1 1 




u s u,yj 
n 


3.2 


where 


S 

u c x >y) -xv 


3.3 
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, A 

For the anti -symmetric part (Fig. 6b) the modes U are chosen as 
follows: consider the classical modal expansion for the pressure 
coefficient for lifting problems 



m-=.o n=o 


which satisfies the Kutta condition at the trailing edge as well as the 
general behavior shown in Fig. 7. By integration with respect to x, 
one obtains (Appendix D) a modal expansion for (5 given by 




3.5 


where $ mn is related to by the relation 


K-= $ 3.6 
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x , n 

(with j = 0) and \J^ n are the desired inodes for the expansion 


for • The expansion for U^is given by 


u «» i*,v> = ± [-JF ci-xj' 


,m + i 


^m j0 w 1 4?i J-4 i -y a v %n 


3.7 


where the plus (minus) sign holds for the upper (lower) surface of 
the wing and 


®tn,o ~ ^ for Vn =. o 
£«t\,o - c T or w ^ o 


3.8 


is the well known Kronecker Delta. 

Combining Eqs. 3.1, 3.2, and 3.5 one obtains 
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9 Q 


Q = U_ (x.y) 

^ i 1 f <i> iYl v 11 


+ U^x, V) 

^ ^ m n tfin 


3.9 


Equation 3.9 may be rewritten as 


<P = X$ A U/WiY) 

A 


3.10 


where 



= 1 


• I 


Pq 


Hi (x,v) - U^WiY) 



,.,PQ 


3.11 
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I 


with 


1 - -f 4- 13 


3.12 


$ . - $ K x-.PQ + 1 , , 

•*“A X m Y\ * * * * 

. .jlPQ+ Kn) 

3.13 

u,. Cx,y) r (x,y) A - Po.'4 -l, . . , 

, ( PO + tANj 


with 


>T — -V- yy> ~ ij 


3.14 



3.3 Formulation 

Consider the equation for the potential on the body surface as 
given by Eq. 1.16 and repeated here for convenience 



3.15 


In the modal method, the unknown is assumed to be a linear 
combination of M known functions, (modes) with M unknown 

coefficients $ Mathematically, this is given by Eq. 3.10 or 

9=1$; W 


Substitution of Eq. 3.16 into Eq. 3.15 yields 
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B 


r 

1 


w 


3.17 


with &^>/c>b given by the boundary conditions. Satisfying Eq. 3.17 at 
M control points yields 


b K + Z E k x^ 
+■ Z W Ka 


where 
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S 


where i 
Equation 3 



3.19 


s defined in Appendix B and is the wake contribution. 
.18 can be rewritten as 


[ A*®] = W 
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3.21 


where the matrix A^.. is given by 

[a,J=M- [ej-M 

Once the unknown coefficients ^ . have been determined, q can be 
evaluated from Eq. 3.16. The coefficient b^can now be evaluated as 

/c)n is known from the boundary conditions. Note tha* once the values 
of £ . are known, the lifting component of the pressure coefficient can 
be evaluated using Eq. 3.4 with K mn given by Eq. 3.6. 

3.4 Least-Square Modal Method 

In the form thus far developed, the modal method can be used but 
is limited to a small number of collocation points, since for more than 
a few collocation points computer times are very large (18 minutes is 
required for 4 collocation points). However, four collocation points 
does not guarantee that the solution is close to the correct one. For 
this reason, a method which allows few modes (4 was found to be sufficient) 
and many collocation points (36 per surface in 18 minutes of computer 
time) is developed here. 

To begin, consider Eq. 3.20 in the following form 
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zAj<£j ~ bx -° 


3.22 


If the number of collocation points, or equations (i) is equal to the 
number of unknowns of modes (j), then Eq. 3.22 can be satisfied 
exactly. Howeve**, for the case when i> j, this equation cannot be 
satisfied exactly, but, rather in the least-square sense. Using the 
method of least-squares, Eq. 3.22 is replaced by 




where E is the sum of the minimum errors. Note here that h 

points 

t Codes' Next » in order to minimize the error E, set the derivatives 
of the error with respect to equal to zero: 
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Using Eq. 3.26 and expanding Eq. 3.25 yields 




Since [A^] = [Aj.] T and rewriting Eq. 3.27 gives 



Finally, for convenience, write Eq. 3.28 as 


[ Q *j]KS 

where 

= [\i] [Aij] 

and 

SbJ=[a«] t W 
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4. NUMERICAL EVALUATION OF THE COEFFICIENTS 


4.1 Introduction 

The purely analytic method. Ref. 3, is limited in that its reliability 
is not assessed for wings of thickness ratios in excess of 1%. This 
is true because of the nature of the method but, on the other hand, a 
purely numerical approach seems well su’'..ed to the problem of thick 
:/’ngs. The final results of the purely numerical method shows that 
although well suited for thick wings, the thicknesses that can effective- 
ly be evaluated are too great to be practical. Considering the results 
of both methods, the reasonable solution is a combination of the two 
methods, which is considered in Section 5. 

4.2 Theory 

First consider the surface ]T divided into elements d£^ as 
shown in Fig. 3. For this case, Eq. 1.15 is rewritten as 





A 


where I is the contribution of the wake as discussed in Appendix B. 
w 

The mean value theorem gives 
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4.2 



where <£ . is the appropriate value in the element In other 
words , Q ^ can be approximated by the value of <?. at the center of 
the box. With this in mind, Eq. 4.1 can be written as 




-£ 


— 1 A T 





X 


4.3 


The 1 in the last term covers only the trailing edge for the contribu- 
tion of the wake with 




4.4 
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where the plus (minus) sign is used for the upper (lower) side of the 

wake and J is as defined in Appendix B. Finally, by satisfying this 
w 

equation at P(k), the center of the box, one obtains 


— b K +>L c ka ^a 4- 


where 



T 


The value of r^ is the distance from the dummy point of integration to 
the canter of the element or 


r K = (0 




JpCvO -p 


Finally W^,. is given by 


4.8 







4.9 


for the boxes along the trailing edge. Boxes other than on the trailing 
edge for wake contribution are given by 



4.10 


The final form for the solution of this problem is given by 
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4.11 


A 


with 





4.12 


where & ^ is the usual Kronecker Delta, viz. 


&\U = O K 4 i 

4.3 Numerical Solution of the Potential Equation 

The potential integral, as given by Eq. 4.1, contains a singularity 
in the ratio 1/r. As the point being evaluated, x, y, z, is approached 
by the dummy point of integration, x-j, y-j, z-j , r tends to zero and the 
integral becomes singular. This problem can be eliminated by the intro- 
duction of polar coordinates 
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dxdy =• j* dr d a 


with the Jacobian given by 

du,y) 

J <JCr,e) 4- 

Also to be considered are the singularities formed by the leading 
edge ( = oo ) and the tip { -y = oo ). These singularities are 

eliminated using the transformation given by Eq. 2.14. A method for 
eliminating both types of singularities simultaneously is shown in 
Appendix C. 

Finally, Eq. 4.6 must be put into numerical form and the Gaussian 
parameters introduced. First consider the source coefficient, b^ 

=JJih H «,■ 
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with the upper surface given by the following: 


S= i -z u U;V) =o 


Rewriting Eq. 4.14 yields 



^LdS 

I VS | 


Introducing Eq. 4.15 gives 


b 




H i ‘ dx t ci Vi 
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4.15 


4.16 


4.17 


r 


i 


Using Eqs. 2.13 and 2.14, Eq. 


★ 

4.17 yields 



'2,'tt r dx 


c 4 X(i _y)dyd y 


For convenience, rewrite Eq. 4.18 as 

b K =jJ flu,/) 


and finally, in suitable numerical form as 


b K HI X x j Yj) 

V — 

Note that the Jacobian of the combined transformation is given by 
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4.18 


4.19 


4.20 




where w-jj, Xj and Tj are the Gaussian quadrature weights and abscissas. 
In a similar manner, consider the doublet coefficients 



4.21 


which can be rewritten as 


C 




S- vr j 


4.22 


Since using Eq. 4.15, iz = dx-j dy^ , the above equation can be rewritten 
as 


c 


*J 




4.23 


using Eqs. 2.13 and 2.14, one has 
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C Kj 


V 5 - v 7 r b. 
~ r* ^ a 


4X (1 ~y )]cxdy 


which, for convenience, can be written as 


4.24 


C * - 

^ KJ 


H (x,y ) dx dy 


4.25 


Finally, Eq. 4.25 can be approximated in numerical form as 


c Kj w ,', H 

J ij 


4.26 
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Finally, t.u system given by Eq. 4.11 can be evaluated to form a 
system, of N linear, simultaneous equations in N unknowns. For a 
high speed digital computer, the generation of this system is limited 

In the number n hv the available comnuter stnraae. The solution can 
be obtained using Gaussian elimination with pivoting. Gaussian 
elimination is quite suitable to the solution of this system as there 
is the possibility of the matrix becoming sinyular. The subroutine 
used to apply the Gaussian method is not particular to the program 
devised but is the IBM subroutine "GEIG" . 

Finally, the solution of the pressure distribution is obtained 
by taking the derivative of (g> with respect to x as shown in Eq. 1.11. 


5. NUMERICAL-ANALYTIC EVALUATION OF THE COEFFICIENTS 


5.1 Introduction 

The purely analytic method gives good results but only for thin 
wings. On the other hand, the purely numerical method gives good 
results but only for thick wings. A good method for evaluating 
moderately thick wings (i.e., wings of thicknesses in the range of 
those in practical use) should incorporate both methods. Tnis 
combination is outlined in this Section. 

5.2 Theory 

This method is based on evaluating the surface both numerically 
and analytically to get the proper results. To begin, the body is 
approximated by tangent planes. To do this, first consider a tangent 
plane element in x, y, z so that z Q is given by 


z - sy 

with the surface, S, given by 

S = 2-«x-?y -2 c = o 
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The exact expression for the doublet contribution is given by 



dy-j 


5.3 


The expression for the doublet contribution in the tangent plane method 
is given by a similar expression as 



V'S-r - vr 



T 



5.4 


with the subscript T indicating a "tangent plant" evaluation* in 
particular 
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The combination of these integrals, the difference integral, becomes 


aC 





(C 





vs > vr 
v 3 



5,6 


This integral is evaluated numerically using the method indicated in 
Section 4. In a similar manner, the correction to the source contribu 
tion is considered as 



2?rr 


aiu 

dX 


_L f iZuj ' 

< 27r'r T ' ax / t 


dx 1 dy 1 


5.7 


also to be evaluated numerically. These differences will be added to 
the tangent plane approximation to yield a better estimate of the 
coefficients Cy and b^. It may be noted that higher numerical errors 
can now be tolerated since "errors on the correction" have small effect 
on the coefficients. 
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6. RESULTS AND CONCLUSIONS 


6.1 Introduction 

The three methods proposed, the modal method, the numerical 

method, and the analytical -numerical method, were proqrammed on a digital 

4 

computer. Comparisons are made against experimental and analytical- 
results. The body surface used is a symmetrical wing at a 5° angle r ' 
attack (see Subsection 2.1). 

6.2 Modal Method Results 

The modal method yields good results but is limited in application 
with respect to the complexity of the body geometry involved. Numerical 
results obtained by using this method were very satisfactory for a good 
range of body thickness ratios. Comparison of the modal method is made 
with the analytical method using the same problem considered in Ref. 3. 

The problem is a rectangular wing in steady subsonic flow (Subsection 2.1). 
The number of modes used for the solution is four. 

First, consider comparison with the box method^. Figures q and 
10 show the comparison between the two methods in the chordwise and 
spanwise directions, respectively. Note that the slope of the diagrams 
for the box method and modal method compare favorably. Rememher that 
the pressure distribution is the derivative of the potential so that 
the importance of the slope agreement can be brought into perspective. 

Next, consider the convergence of the medal method in terms of the 
number of evaluation points. Figures 11 and 12 depict the convergence 
in terms of the number of points of evaluation on the body surface in the 


- 47 - 



chordwise and spanwise directions, respectively. Agreement here is 
extremely good. It may be noted that the solution using this metnod is a 
smooth curve and not discrete points. Points between the solution 
for six points show excellent agreement for five, four and even three 
points. Also, it is observed that fewer points describe the curve 
with less accuracy, but even for few points the agreement is still 
reasonably good. 

Finally, the effect of the thickness ratio is considered. In 
Figs. 13 and 14 the convergence for hody thickness ratios of .1, .01, 
and .001 is shown in the chordwise and spanwise directions, respective- 
ly. The results here are typical of those for the box method^. Again, 
the behavior here is what is expected as discussed in Subsection 6.2 

Results for the modal method show it to be a reliable method 
for varying thickness ratios but, unfortunately limited in its use due 
to the difficulty in “guessing" the shapes of the potential, and thus 
the modes, for complex geometries (i.e., wing-body interference). 
Furthermore, the computer time is much longer than for the box method. 
Hence, as mentioned above, the box method is considered to be superior 
and was investigated more deeply. 

6.3 Analytic-Numerical Results 

The analytical-numerical method provides very satisfactory results 
in a varying range of thickness ratios. It is the varying range of 
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applicability with respect to body thicknesses that makes this method 
sc desirable. 

The c* ^ a .yticc '-numerical method is compared here with experiments 
tnd also with the purely analytical method in order to evaluate 
the improvement. Furthermore, the convergence of the method and the 
thickness effects are also examined. 

Figure 15 shows a comparison of the analytical method, the analytical 
numerical method and experiments using the wing in Subsection 2.1 at 

O 

an angle of attack of 5 and a thickness ratio of 1%. Results here show 
that the analytical -numerical method provides a good description of the 
pressure distribution with respect to exDeriments and therefore, the 
lift, the area under the pressure diaqram, is also described well. Note 
also that both the analytical and the analytical -numerical methods compare 
favorably. 

Next consider the convergence with resDect to the numher of points 
of solution on the body surface. Since the solution is a solution of 
discrete points it is necessary to check that when the number of 
points of evaluation is changed, the solution converges to the same 
curve. This is shown in Figs. 16 and 17 for various numbers of points 
in the spanwise and chordwise directions, respectively. For a wing with 
six points of evaluation connected by a smooth curve it is evident 
from Figs. 16 and 17 that the solution based on 3, 4 and 5 points lies 
satisfactorily on the same curve. 
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Finally, consider the influence of body thickness on the results, 
Here the potential is plotted in tue spanwise and chorawise directions 
in Figs. 18 *nd 19, respectively. The thickness ratio is varied in 
these figures fr ;m .001 to .1. The behavior is as expected. Note 
that at the tip and leading edge, the variation is *mall, while at the 
root and trailing edge the variation is larger. 

In summary, it is evident that the analytical -numerical method 
provides good results for a large range of thickness ratios. Note 
that a body thickness ratio of .1 is a wing with maximum body thickness 
to wing chord ratio cf 10%, well within practical body thickness ratios. 
Since the method is also convergent, it has proven to be quite a re- 
liable and accurate method. 

6.4 Suggestions for Further Study 

The geometry evaluated in this thesis are of the most simple 
case, a symmetrical wing at a small angle of attack. Obviously, this 
geometry was chosen for its simplicity in testing the methods as it 
is easier to program and debug. It is also quite obvious that the most 
advantageous quality of this method is the evaluation of complicated 
geometries to include the problem of wing-body interference. It is 
therefore suggested that a study of complicated geometries be made using 
the analytical -numerical mettod. 
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Because of the nature of the formulation of the problem, it is by 
no means restricted to steady. Incompressible flow as considered by this 
thesis. It is suggested that the numerical formulation be extended to 
include unsteady compressible (subsonic and supersonic) flow-fields. 
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APPENDIX A 


DERIVATION OF THE INTEGRAL POTENTIAL EQUATION ON THE BODY SURFACE 

For the sake of completeness in this appendix, the derivation 
of Eq. 1.12 (Ref. 3) is summarized. For convenience, Eq. 1.12 is 
given here as: 




^ ati/ x, ' Y i > z ^4Trr 





A.l 


where 


r 





A. 2 


When the point P is allowed to approach the point P*, as shown in 
Fig. A1 , the integrands of Eq. A.l become singular. In order to per- 
form the limit it is convenient to isolate a neighborhood about P* 
having area cr* and radius £ . With this in mind, rewrite Eq. A.l 

as: 
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r r ^<r -» . — 

j>i j^) — ~ ^js cJX 

x-<r 

I + jf 

cr x-<r 

+ - <>..(■,.?,) |- (WFJ dI 


- 1; - 1; < *i*. 


H ii 


With obvious definitions of I s , 1$ and I Q . 
First consider the doublet integral 

r = i' + i" 

D t> S 


A. 3 


Consider a local coordinate system ^ , 5 as shown in Fig. A2 

in order to evaluate I Q ". Assume the area cr to be circular having 
radius £ , so small that <r can be replaced by its tangent plane. 

tl 

Rewriting I Q in the new variables yields: 
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It 

O O 


V° 


3- _ 

where £ = + n, and o< ii shown in Fig. A3. Note that 

r a = (X, - x)"+ (y, -y/ Kz^Ej a 


* a 


= s i + n , +( s t * s *v 

= R a +[(5 1 -*)* 


v° 


v° 


Next consider 


^ f _L_ * 

<ss, v 4'irt' j 


i 

lv 


1- ( 5 ' ~ ^ 
2 ^ TT r 3 


i t,-o 


4 ?rH 


and 


4, ( — ! \ = _ J_ a 0 

,J 9 \ 4?rr j a 4^3 

= - 9 

4Tf 


Combining yields 






'O 


— ^ f 1 


Using A. 6, Eq. A. 4 becomes 


// 






C o 


a. 
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For very small £ , <p is nearly constant and can be replaced under 
this assumption by its value at the point x*, y*, z* as 


9*- 


Equation A. 7 now becomes 


I " = - S 


O 






- Airr J 




»- o 


A. 8 


) 


= ( i - - 5 — 

As y goes to zero, i.e., x, y, z approaches x*, y*, z*, Eq. A. 8 
becomes 


livn 


-i n 




•Iffi 
5 50 
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; i - 


-1 £ a +r, 
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A. 9 


Combining Eqs. A. 3 and A. 9 and applying the limit as P-* P* yields 




UT = «* 

z 


+" li m ' <f£ (ft ( —1— Ur 

r-<r c 


A. 10 


with 


t* = [ (** - X. 4- (y K - ^f+- (z x -z/j 


y* 
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Similar treatment of the sources yields no special contribution from 
the neighborhood of P*. Equation A.l can now be written as: 



After combining terms, the form of Eq. A.l on the body surface is given 
as 



1 I y 


-f* 



A. n 


where the second integral is obtained as a limiting case for£~>0. 
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APPENDIX B 


WAKE FORMULATION 


For completeness, the formulation of the wake given in Ref. 3 is 
summarized here. According to the linearized Bernoulli Equation, 

Eq. 1.15, for steady flow, there can be no pressure jump through the 
wake so that on the wake 


o>X< 




= O 


B.l 


For Eq. 1.15, the expression for I w is written as 



B.2 


where 2T w is the upper surface of the wake and n-j is understood to be 
the upper normal. The surface 7L w is shown in Fig. B1 . Note that the 
surface of the wake is not known and Eq. 1.14, which is satisfied on the 
body, must be completed by the equation on the wake. From Eq. 1.15, 
the velocity or. the wake is tangent to the surface of the wake. What 
results is two coupled equations, one on the body with c)Q/c>tn known 
and ^ unknown and the second on the wake with /dr\ unknown and 

( =< ? A - (Q i ) constant along the x - direction on the wake and 


equal to the value of A<p at the trailing edge. The geometry of the 
wake is determined given the condition that the velocity is tangent to 
the wake by knowing the velocity on the wake, d<?/dn . 

To simplify the above procedure, consider the following. First 
note that 



where 


dF„ = = dF.c.os(t>,r) 


and is the solid ergle shown in Fig. 26 as 
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Next, consider the wake integral as a sum of strips of width dy in the 
x direction as shown in Fig. B2. Applying the mean value theorem to 
Eq. B.2 combined with B.3 yields 


Iw = 


- £<F(y; d-fl- 

<- f( 

- 2.AQKVJ) ! : JjV 

j Ju 


B.5 


By noting that the Integral Is the solid angle for each element 
^•j one obtains 


T = j B.6 

i VV J 

(Note that a<£ is a function of y only) 

From Eq. B.6, it can be seen that the wake contribution is u 
function of the solid angle so that a reasonable geometry for the 
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wake is one in which the solid angles remain fairly constant. This 
suggests the possibility of approximating the wake by straight vertex 
lines, parallel to the direction of flow, emanating from the trailing 
edge of the wing since geometrical considerations show that the solid 
angles change only slightly? With this assumption, the wake equation 
simplifies considerably and its contribution reduces to a line integral. 
Assuming the trailing edge to be given by 

x = 

2 - 2 . rt (y) 

the equation of the wake surface becomes 

S w r Z, -Z-rJVj 


and 


I w = -jjtSSf ( VS • r)-p 

- V) 4- Zvi (Vi) 

~/ A<? ^ d Vl 
-% 


\ . 
i 

V 

y 
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where b is the span of the wing and 



with Z.j * Z W (Y). If the trailing edge is in the plane Z«j= 0, 
Eq. B.8 reduces to 

T = — ±JL 1 i - 

w W^V) 2 + Z 2 _ 

X te ~ X 

+W,-Y) 2 *• H 2 y iz 


Note that the above discussion is for bodies with a sharp trailing 
edge only as it depends upon the known location of the stagnation point 
from which the wake emanates. For bodies without a sharp edge, e.g.» 
a rotating cylinder, this formulation is not applicable. 
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APPENDIX C 


As mentioned in Subsection 5.2, the coefficients 

-Jfk'' FSW 

r x ' 

are evaluated by approximating the surface element with its tangent 

plane at the centroid of the element. This value can eventually be 
corrected by adding the integral of the difference between the original 
integrand and the (tangent-element) approximate integrand. This 
integral can be evaluated numerically by using standard Gaussian 
quadrature formulas. However, in the case of k = 1 (effect of the 
element on itself) the tangent plane contribution is equal to zero. 
Furthermore, the integrand becomes infinite when = 0. Hence, a 
special integration scheme must be used. In this appendix an analysis 
of the type of singularity of the integrand of Eq. C.l (when k = i) 
is given. Then a transformation that eliminates the singularity is 
presented. 

First consider the behavior of the doublets in the neighbornood 
of the singularity. For simplicity, the analysis of the behavior of the 
doublet in the neighborhood R s 0 is performed with a frame of reference 
such that origin is at the centroid of the element and the z axis is 
directed as the normal n. Then the equation of the element can be 


y 
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written as 


Z - FU.y^ =o 


C.2 


with z = 0forx = y = 0 (Fig. C.l) and Eq. C.l for k = i reduces 
to 



C.3 


where 


h = " a* lx i' ^ ” Hi (y '~ y } + (z; ‘ z) 


C.4 


is the distance along the normal of the origin from the point Xj , 

, z^, (Fig. C.2). If R goes to zero the distance h goes also to zero. 
More precisely, as R goes to zero 


ORIGINAL PAGE IS 
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Vi Q R z c 

where Q = 1/^- is the curvature of the cross section indicated in the 
figure is the radius of curvature). Thus, in a neighborhood of 

R * 0, Eq. C.3 reduces to 


c 





dy ( 


C.6 


It should be noted that Q is the curvature of the cross section and 
thus Q depends upon the angle of the cross section. Thus, in 

order to evaluate Eq. C.3, it is convenient to use polar coordinates 
since this eliminates the singularity as well as the sharp variations 
(in the plane, x^ » y^) of the integrand due to the dependence of Q 
upon ^ 

Next consider the integration scheme. As shown above, the 
Integration of Eq. C.l (with k = i) In the neighborhood of the 
centroid of the element can be performed using standard quadrature 
techniques (Gaussian quadrature in particular) in polar coordinates. 
However, the domain of integration is not simply defined with these 
variables. Hence, a more suitable technique (fully correspondent 
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to Integration in polar coordinate) is described here. 

As mentioned in Subsection 2.3, the boundary of the element in 
the plane are given by 

x: - a* £ ^ £ x c + ax 

c 

C.7 

y c - Ay £ y £ y c + Ay 


Note that the use of x, y variables 

. -1 
^ = X 

n =[i-u-v) 1 ] ^ 


has the advantage of eliminating the square root singularity at the 
leading edge and the tip*. On the other hand, a singularity of the 
type K ^ in the plane x-j , y^ yields a singularity of the type R 
in the plane x, y. Thus the integral to be evaluated is of the type 



*Th1s singularity is due to the factor 17 s | in d£ c. 
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where f(x, y) Is a finite but discontinuous function of x and y (the 
discontinuity being due to the “cross-section-curvature effect") but 
continuous in polar coordinates. 

In order to analyze Eq. C.7, it is convenient to separate the 
contributions of and (Fig. C.3) as 


with 



c.n 


c kk = jj k x < y) dxcl >' 

1* t>4 

Using the transformation 

X = X t +- ~ u -1 £ U 4 +1 

y - Y + ITU ~ < - V - tl 

c / 


C.12 


C.13 
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Equation C.ll reduces to 


+i 

= 

Note that f(u, v,)u is a regular function of u and v since the factor 
u compensates for the (x^ + y2) ^ singularity and the "cross- 
section-curvature effect" disappears in the u, v plane (which is similar 
to polar coordinate plane). Hence, Eq. C.14 can be evaluated by Gaussian 
quadrature. A similar transformation can be used in Eq. C.12. 

This procedure was used to evaluate not only the effect c^ of 
the element on itself but also for the effect of an element on the 
opposite element (e.g., leading edge root boxes on opposite surfaces). 

A similar technique is also used for the evaluation of 



k Ki 
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APPENDIX D 


MODAL METHOD DERIVATIONS 


This Appendix gives the derivations of the mode shapes for the 
symmetrical and anti -symmetrical modes and the expression for the 
pressure coefficient. First, consider the symmetrical modes as given by 

W) Oj 


From Fig. 6A, the symmetric components of the potential can be approxima- 
ted in the chordwise direction by 

U«00 = X n *=0,1, D.2 

and the spanwise direction by 

\J m { y> =y 'ff'=o,i v z , . . . o.3 

Next to be considered is the anti -symmetrical modes as shown in 
Fig. 6B. Here not only the shape of the potential, but also the shape 
of its derivative which has the form of the well known pressure distri- 
bution curve must also be derived. Note that the curve is infinite at 
the leading edge and drops to zero as the square root of x, the distance 
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from the trailing edge. Consequently, a function satisfying all these 
conditions is necessary. A function such as this is classically given 
by 

C p =-2<P* 




D.4 


•'here represents the derivative of with respect to x. Examination 

of Eq. D.4 shews that it gives the required conditions at the leading 

edge, where x = 0, and is singular (i.e., C - OO ) and also that as x 

P 

approaches 1, near the trailing edge, it behaves as the square root 
of the distance from the trailing edge, 1-x. Since Eq. D.4 satisfies 
all the necessary conditions, it can be integrated to obtain <£ or 



where the limits are from the leading edge to the trailing edge, c. 
Integration of this equation is shown in the following: 
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where n ~ m + 1. Combining the first and third terms yields 

1= ^^(^(1-x) 

< 1 - xf] - 

Next, integrating 
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Finally, consider the formulation for the pressure distribution. 

From Eq. 1.11, the expression relating C with the potential involves 

P 

a derivative with respect to x, the direction of the flow. From this 
relation, it is seen that the expression for the modes is an analytical 
differentiation with respect to x. 

Begin by considering the symmetric component of the potential 
in the chordwise direction only. From Section 3 this is given by 

U* (X) - x" 

The pressure coefficient mode shape resulting from this component 
is given by 


a 

ax 
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Next, the anti -symmetric component is considered. Care must be 
takpr as the arc -tangent term is considered only n = 0, the first mode. 
The upper-surface component is given as 
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First, consider the evaluation of pressure coefficient, Eq. 1.11, 

for the kth mode, U, , Eq. 0.12. 
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u' K = fxMi-w*'* 
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+ -fx U-X) VK-tH-1) 


which becomes 
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It should be noted that an arc-tangent term appears in the first mode 
(see Eq. D.12). Differentiating the first mode yields 
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Note that the last term contains the arc-tangent contribution. Combining 
terms gives 


- 98 - 


f 


-i -i ] 

or 


u 


t 

1 


i^ftZ* 

2 V x 


D.15 


Equations D.14 and D.12 give the relations for the pressure coefficient 
for the first through Kth codes. 
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WAKE GEOMETRY 
FIGURE 81 



WAKE EVALUATION 
FIGURE 82 




DOUBLET BEHAVIOR IN THE AREA OF THE SINGULARITY 

FIGURE Cl 
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